Fastfood restaurants by cuisine type

Fastfood restaurants are identified by their cuisine descriptions given in the inspection data. We print out the cuisine descriptions list (n=85) and let everyone circle the ones they think is fastfood and the union are used as our rule (use union because it’s more conservative).

We classify cuisine descriptions “Bagels/Pretzels”, “Bottled beverages, including water, sodas, juices, etc.”, “Chicken”, “Delicatessen”, “Donuts”, “Hamburgers”, “Hotdogs”, “Hotdogs/Pretzels”, “Ice Cream, Gelato, Yogurt, Ices”, “Nuts/Confectionary”, “Pancakes/Waffles”, “Pizza”, “Soul Food”, “Sandwiches”, “Sandwiches/Salads/Mixed Buffet” and “Soups & Sandwiches” as fastfood restaurants. And then we calculate the total number of restaurants and the number of fastfood restaurants, as well as the percentage of fastfood restaurants for each neighborhood.

# calculating the total number of restaurants and the number of fastfood restaurants in the neighborhood, as well as the percentage of fastfood restaurants.
neighborhood_list = 
  rest_neighborhood %>%
  distinct(neighborhood) %>%
  arrange(neighborhood)
  
rest_fastfood_neighborhood = 
  rest_neighborhood %>%
  filter(cuisine_description %in% c("Bagels/Pretzels",
                                    "Bottled beverages, including water, sodas, juices, etc.",
                                    "Chicken",
                                    "Delicatessen",
                                    "Donuts",
                                    "Hamburgers",
                                    "Hotdogs",
                                    "Hotdogs/Pretzels",
                                    "Ice Cream, Gelato, Yogurt, Ices",
                                    "Nuts/Confectionary",
                                    "Pancakes/Waffles",
                                    "Pizza",
                                    "Soul Food",
                                    "Sandwiches",
                                    "Sandwiches/Salads/Mixed Buffet",
                                    "Soups & Sandwiches"))

percent_fastfood_neighborhood = function(name_neighborhood){

  rest_each_neighborhood =
    rest_neighborhood %>%
    filter(neighborhood == name_neighborhood) %>%
    distinct(camis)
  
  n_rest_neighborhood = nrow(rest_each_neighborhood)

  rest_fastfood_distinct_neighborhood = 
    rest_fastfood_neighborhood %>%
    filter(neighborhood == name_neighborhood) %>%
    distinct(camis, cuisine_description)
  
  n_fastfood_neighborhood = nrow(rest_fastfood_distinct_neighborhood)
    
  percent_fastfood_neighborhood = n_fastfood_neighborhood/n_rest_neighborhood
  
  tibble(
    neighborhood = name_neighborhood,
    n_fastfood = n_fastfood_neighborhood,
    n_rest = n_rest_neighborhood,
    percent_fastfood = percent_fastfood_neighborhood
  )
}

fastfood_neighborhood = 
  map(neighborhood_list$neighborhood, percent_fastfood_neighborhood) %>%
  bind_rows() %>%
  mutate(neighborhood = str_to_upper(neighborhood))

# plot for each neighborhood
fastfood_neighborhood %>%
  mutate(neighborhood = as.factor(neighborhood),
         n_rest = as.numeric(n_rest),
n_nonfastfood = n_rest - n_fastfood,
         neighborhood = fct_reorder(neighborhood, percent_fastfood)) %>%
  plot_ly(., x = ~neighborhood, y = ~n_fastfood, type = 'bar', name = 'fastfood restaurants') %>%
  add_trace(y = ~n_nonfastfood, name = 'non-fastfood restaurants') %>%
  layout(yaxis = list(title = 'Number of restaurants'), 
         xaxis = list(title = 'Neighborhood (ordered by percentage of fastfood restaurants)',
                      showticklabels = FALSE),
         barmode = 'stack')

While the Greenwich Village and SOHO neighborhood has fairly large number of restaurants, it has the smallest percentage of fastfood restaurants. Williamsbridge and Baychester has the largest percentage of fastfood restaurants.

When large number of total restaurants is not equal to large percentage of fastfood restaurants in that neighborhood, we can conclude that the distribution of fastfood restaurants is not even across neighborhoods, which also implies the motivation of our study, we want to investigate if this uneven distribution of fastfood restaurants is associated with diffferent level of chronic disease outcomes within a neighborhood.

Restaurant Chains

We first scrape the list of 75 national chain restaurants in the US from the wikipedia page (https://en.wikipedia.org/wiki/List_of_restaurant_chains_in_the_United_States#Fast-casual) and then join this dataset with restaurant inspection data to choose only the chain restaurants in NYC.

chains_html = read_html("https://en.wikipedia.org/wiki/List_of_restaurant_chains_in_the_United_States#Fast-casual")

# read in the list of chain restaurants in us 
# made the names to uppercase and changed the var name to dba
chain_rest = chains_html %>%
  html_nodes("td:nth-child(1)") %>%
  html_text() %>%
  as.tibble() %>%
  mutate(dba = value, 
         dba = str_to_upper(dba),
         dba = str_replace_all(dba, "[\r\n]" , "")) %>%
  select(dba)

## Warning: `as.tibble()` is deprecated, use `as_tibble()` (but mind the new semantics).
## This warning is displayed once per session.

head(chain_rest, 10)

## # A tibble: 10 x 1
##    dba                               
##    <chr>                             
##  1 ""                                
##  2 A&W RESTAURANTS                   
##  3 AMERICA'S INCREDIBLE PIZZA COMPANY
##  4 APPLEBEE'S                        
##  5 ARBY'S                            
##  6 ARCTIC CIRCLE RESTAURANTS         
##  7 ARTHUR TREACHER'S                 
##  8 ATLANTA BREAD COMPANY             
##  9 AUNTIE ANNE'S                     
## 10 BAHAMA BREEZE

Then, we match the list of chain restaurants in U.S. with the restaurant inspection data.

# removing punctuations in chain_rest & restaurant inspections (neighborhoods)
chain_rest_str = 
  chain_rest  %>%
  mutate(dba = str_replace_all(dba, "[[:punct:]]", ""))

rest_neigh_str = rest_neighborhood %>% 
  mutate(dba = str_replace_all(dba, "[[:punct:]]", ""))

# Matching the two datasets(restaurant inspection data that has all punctuation removed from dba(restaurant name) and list of chain restaurants by dba)
neighborhood_chain = 
  right_join(rest_neigh_str, chain_rest_str) %>%
  filter(!is.na(camis)) %>%
  distinct(camis, dba, neighborhood, boro)

## Joining, by = "dba"

neighborhood_chain %>%
  group_by(dba) %>%
  summarise(n = n()) %>%
  arrange(desc(n)) %>%
  head(10)

## # A tibble: 10 x 2
##    dba                        n
##    <chr>                  <int>
##  1 STARBUCKS                286
##  2 SUBWAY                   285
##  3 MCDONALDS                205
##  4 POPEYES                  101
##  5 BURGER KING               78
##  6 CHIPOTLE MEXICAN GRILL    76
##  7 DUNKIN DONUTS             51
##  8 WENDYS                    43
##  9 KFC                       36
## 10 LITTLE CAESARS            35

The combined dataset “neighborhood_chain” has 1546 observations. Also, there were 64 different chain restaurants extracted.

# counting chains in neighborhoods
neigh_count_chain = neighborhood_chain %>%
  group_by(neighborhood, boro) %>%
  summarise(chain_n = n())

neigh_count_rest = rest_neighborhood %>%
  distinct(neighborhood, camis) %>%
  group_by(neighborhood) %>%
  summarise(res_n = n())

# calculating percentage of chains in each neighborhood
percent_neighborhood_chain = left_join(neigh_count_chain, neigh_count_rest) %>%
  ungroup() %>%
  mutate(chain_percentage = chain_n/res_n,
         neighborhood = str_to_upper(neighborhood))

## Joining, by = "neighborhood"

plot_chain_neighbor = percent_neighborhood_chain %>%
  mutate(neighborhood = forcats::fct_reorder(neighborhood, chain_percentage)) %>%
  ggplot(aes(neighborhood, chain_percentage, fill = boro)) + geom_bar(stat="identity") +
  labs(x = "Neighborhoods in NYC", y = "Percentage of chain restaurants") +
  theme(axis.text.x = element_blank(), axis.ticks = element_blank())

ggplotly(plot_chain_neighbor)

max(percent_neighborhood_chain$chain_percentage)

## [1] 0.1748252

We plot neighborhoods in NYC with their percentage of chain restaurants and group them by borough. We can see that neighborhoods with those smallest percentages of chain restuarants are mostly in Brooklyn except for Greenwhich Village and Soho and Lower East Side and Chinatown in Manhattan. Neighborhoods in Queens and Manhattan are spread out across low to high percentage of chain restaurants while most of the neighborhoods in Bronx and Staten Island have high percentages. The neighborhood with the highest percentage of chain restaurants is Throgs neck and Co-op City in Bronx with 17.5% of chain restaurats out of all restaurats.

Inspection Grade

gradea_neighborhood =
  rest_neighborhood %>%
  group_by(neighborhood, grade) %>%
  summarise(n = n()) %>%
  mutate(grade_percent = n / sum(n)) %>% 
  filter(grade == "A") %>% 
  ungroup(boro) %>%
  mutate(neighborhood = str_to_upper(neighborhood))

gradea_neighborhood %>% 
  mutate(neighborhood = as.factor(neighborhood),
         neighborhood = fct_reorder(neighborhood,grade_percent)) %>% 
  plot_ly(x = ~neighborhood, y = ~grade_percent, color = ~neighborhood, type = "bar") %>%
  layout(yaxis = list(title = 'Percentage of Grade A restaurants'), 
xaxis = list(title = 'Neighborhoods in NYC', showticklabels = FALSE),
         showlegend = FALSE) 

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

The differences on percentage of “grade-A” restaurants between each neighborhood are observed. “Throgs Neck and Co-op City” has the greatest grade-A restaurant percentage, around 50.7%. “Sunset Park”, however, has the least, around 32.4%. “Washington Heights”, where we live, takes the fourth counting from the end, around 34%, which is obviously consistent with the feeling we have towards the restaurant condition of “Washington Heights”

health_neighborhood = 
  health %>%
  mutate(neighborhood = str_to_upper(name)) %>%
  select(-name)

combined_chain = 
  percent_neighborhood_chain %>%
  select(neighborhood, chain_percentage)
combined_chain_fastfood = 
  fastfood_neighborhood %>%
  mutate(fastfood_percent = percent_fastfood) %>%
  select(neighborhood, fastfood_percent) %>%
  right_join(combined_chain, by = "neighborhood")
combined_chain_fastfood_gradea = 
  gradea_neighborhood %>%
  select(neighborhood, grade_percent) %>%
  right_join(combined_chain_fastfood, by = "neighborhood")

combined_model =  
  left_join(combined_chain_fastfood_gradea, health_neighborhood, by = "neighborhood")

Main predictor selection

outcome_name = combined_model[,10:12]

main_predictor_selection = function(outcome){

lm_fastfood =   
  lm(outcome ~ fastfood_percent + grade_percent + racewhite_rate + poverty + smoking + exercise + airquality_rate, combined_model) %>%
  summary() 
lm_fastfood_tbl = 
  as.tibble(lm_fastfood[[4]]) %>%
  clean_names()

lm_chain =   
  lm(outcome ~ chain_percentage + grade_percent + racewhite_rate + poverty + smoking + exercise + airquality_rate, combined_model) %>%
  summary() 
lm_chain_tbl = 
  as.tibble(lm_chain[[4]]) %>%
  clean_names()

lm_both =   
  lm(outcome ~ fastfood_percent + chain_percentage + grade_percent + racewhite_rate + poverty + smoking + exercise + airquality_rate, combined_model) %>%
  summary() 
lm_both_tbl = 
  as.tibble(lm_both[[4]]) %>%
  clean_names()

tibble(p_fastfood_sing = lm_fastfood_tbl$pr_t[2],
       p_chain_sing = lm_chain_tbl$pr_t[2],
       p_fastfood_both = lm_both_tbl$pr_t[2],
       p_chain_both = lm_both_tbl$pr_t[3])
}

main_predictor_comp = 
  map(outcome_name, main_predictor_selection) %>%
  bind_rows()
main_predictor_comp$outcome = colnames(combined_model)[10:12]
main_predictor_comp

## # A tibble: 3 x 5
##   p_fastfood_sing p_chain_sing p_fastfood_both p_chain_both outcome    
##             <dbl>        <dbl>           <dbl>        <dbl> <chr>      
## 1        0.000305      0.00221          0.0370        0.387 obesity    
## 2        0.0313        0.0609           0.239         0.594 diabetes   
## 3        0.0137        0.340            0.0126        0.280 stroke_hosp

We originally have two main predictors of interest, Percentage of chain restaurants and Percentage of fastfood restaurants, and they are both significantly associated with the three chronic disease health oucomes when solely in the model after adjusting for other potential confounders. Taking obesity as example, the p-value for fastfood_percent is 0.00030506 and for chain_percentage is 0.0022107.

And we also anticipated them to be highly associated with each other, to avoid collinearity, we need to make decision on which one to keep as the final main predictor. So we put the two main predictor candidates in the same model and see how their p-values change. As a result, the fastfood_percent stays significant (p=0.0370412) and the chain_percentage turn insignificant (p=0.3867866).

Same selecting processes are repeated for the oucomes diabetes and stroke. Thus, we decide to have Percentage of fastfood restaurants (fastfood_percent) as the main predictor.

Confounder selection

outcome_name = combined_model[,10:12]

confounder_percent_change = function(outcome){
lm_adjusted = 
  lm(outcome ~ fastfood_percent + grade_percent + racewhite_rate + poverty + smoking + exercise + airquality_rate, combined_model) %>%
  summary() 
lm_adjusted_tbl = 
  as.tibble(lm_adjusted[[4]]) %>%
  clean_names()

lm_crude =
  lm(outcome ~ fastfood_percent + racewhite_rate + poverty + smoking + exercise + airquality_rate, combined_model) %>%
  summary() 
lm_crude_tbl = 
  as.tibble(lm_crude[[4]]) %>%
  clean_names()

percent_change = (lm_crude_tbl$estimate[2] - lm_adjusted_tbl$estimate[2]) /lm_crude_tbl$estimate[2]

confounder_percent_change = percent_change
}

confounder_change =
  map(outcome_name, confounder_percent_change) %>%
  bind_rows()

percent <- function(x) {
  paste0(formatC(100 * x, digits = 3), "%")
}

confounder_change %>%
  kable()

obesity	diabetes	stroke_hosp
-0.0921697	0.2017842	-0.111628

Besides some biologically meaningful covariates (i.e. race, poverty, smoking status, exercise), we also hypothesize the variable Percentage of grade A restuarants as a potential confounder in the association between fastfood restaurants percentage and the three chronic disease health outcomes.

Here, we are assessing if grade_percent is a significant confounder. Using the 10% change rule of thumb, we find that after adjusting for Percentage of grade A restuarants, the estimates of fastfood_percent change by -9.22% for outcome obsesity, 20.2% for outcome diabetes and -11.2% for outcome stroke hospitalization.

Thus, regarding the final models, we are going to keep grade_percent for obsesity and stroke hospitalization, but take it out and rerun the model for the outcome diabetes.